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Abstract 

The present paper deals with the numerical study of high pressure H2/L02 combustion for propulsion 
systems. The present research effort is driven by the continued interest in achieving low cost, reliable access to 
space and more recently, by the renewed interest in hypersonic transportation systems capable of reducing 
time-to-destination. Moreover, combustion at high pressure has been assumed as a key issue to achieve 
better propulsive performance and lower environmental impact, as long as the replacement of hydrogen with 
a hydrocarbon, to reduce the costs related to ground operations (propellant handling, infrastructure and 
procedures) and increase flexibility. 

Starting from this background, the current work provides a model for the numerical simulation of high- 
pressure turbulent combustion employing detailed chemistry description, embedded in a Reynolds averaged 
Navier-Stokes equations solver with a Low Reynolds number k — to turbulence model. The model used 
to study such a combustion phenomenon is an extension of the standard flamelet-progress-variable (FPV) 
turbulent combustion model combined with a Reynolds Averaged Navier-Stokes equation Solver (RANS). 
In the FPV model, all of the thermo-chemical quantities are evaluated by evolving the mixture fraction Z 
and a progress variable C. When using a turbulence model in conjunction with FPV model, a probability 
density function (PDF) is required to evaluate statistical averages (e.g., Favre average) of chemical quantities. 
The choice of such PDF must be a compromise between computational costs and accuracy level. State- 
of-the-art FPV models are built presuming the functional shape of the joint PDF of Z and C in order to 
evaluate Favre-averages of thermodynamic quantities. The model here proposed evaluates the most probable 
joint distribution of Z and C without any assumption on their behavior with the Statistically Most Likely 
Distribution (SMLD) framework. This provides a more general model in the context of FPV approach. 

I. Introduction 

The development of new technologies to enhance and control combustion processes is nowadays fun¬ 
damental in order to increase efficiency and reduce emissions in many engineering applications, such as 

*Ph.D Student, a.coclite@poliba.it 

'I'Research engineer, Propulsion Unit, AIAA Senior Member, l.cutrone@cira.it, +390823623108 
^Professor, AIAA member, depalma@poliba.it, +390805963226 
^Professor, pascazio@poliba.it, +390805963221 


1 of[T6] 


American Institute of Aeronautics and Astronautics 




reciprocating internal combustion engines; advanced gas turbine systems, pulse detonation engines, high¬ 
speed air-breathing propulsion devices. Hydrogen is one of the preferred fuel because of its properties in 
terms of very short ignition delay time and high energy per unit weight. The investigation of hydrogen 
combustion presents significant difficulties and high costs either following the experimental approach or the 
numerical one. Moreover, high-Reynolds-number turbulent combustion is a formidable multi-scale problem, 
where the interaction between chemical kinetics, molecular, and turbulent transport occurs over a wide range 
of length and time scales. These features pose severe difficulties in the analysis and comprehension of the 
basic phenomena involved in supersonic combustion. 

In this context, the simulation of turbulent reacting flows is very useful to cut down experimental costs 
and to advance the comprehension of the basic physical mechanisms. Turbulent combustion is a multi-scale 
problem, where the interaction between chemical kinetics, molecular, and turbulent transport occurs over 
a wide range of length and time scales. The numerical simulation of such phenomena with detailed chem¬ 
istry is today still prohibitive. The huge computational cost stems from the fact that, even for a simple 
fuel, detailed kinetic mechanism can involve thousands of reactions and, consequently, hundreds of chemical 
species. In recent years, the need for efficient tools has driven the research towards: i) developing models 
for turbulent combustion in order to understand and accurately mimic the interaction between turbulence 
and chemistry ii) studying improved kinetic schemes to describe the combustion process [5HH]- More¬ 
over, simplified approaches to combustion modeling have been proposed to further reduce the number of 
equations to be solved; for instance, the reduction of the chemical scheme in intrinsic low dimensional man¬ 
ifolds (ILDM) [^; the flamelet-based approaches such as the flamelet-progress-variable (FPV) [TU] or flame 
prolongation of ILDM (FPI) [TT] ; and Flamelet Generated Manifolds approach (FGM) [T^] . 

An additional modeling difficulty is related to real-gas effects in high-pressure conditions. For high- 
pressure propulsion systems, the pressure of the injected fluid is always supercritical, whereas its temperature 
may be either subcritical or supercritical. However, the fluid is typically injected into a chamber in which 
both the pressure and the temperature exceed the critical values for both the fuel and the oxidizer, so 
that a fast transition to the supercritical state is observed m, and the liquid phase can be described as a 
dense gaseous jet. At such operating conditions, the ideal-gas equation of state cannot predict the correct 
p-v-t relation for the oxidizer and the fuel; for example, the density of oxygen predicted using the ideal- 
gas equation of state at supercritical conditions can be one fourth of the real value. Moreover, real-gas 
effects have a significant impact on the flame structure in high-pressure combustion, as shown in Ref m 
for the case of a premixed hydrogen/oxygen flame, and in Ref (IlSp for that of counterflow hydrogen/oxygen 
diffusion flame. Therefore, a suitable equation of state, together with adequate constitutive equations for 
the transport properties, may be warranted. 

The present work is based on the FPV model for non-premixed turbulent flames with a particular choice 
for the probability density function (PDF) needed to evaluate the thermo-chemical Favre averages combined 
with an an Eulerian single-phase numerical method for computing the mixing and combustion of liquid 
propellants at operating conditions typical of rocket combustion chambers which computes the detailed 
chemistry for cryogenic applications efficiently and introduces a real-gas model in the flamelet library, so as 
to compute high-pressure combustion accurately m- The aim of this work is to study the applicability of 
the statistically most likely distribution (SMLD) [T7] approach to model joint-PDF of Z and A in the case of 
high-pressure combustion. The proposed joint-SMLD approach is very interesting since it represents a good 
compromise between computational costs and accuracy level. The results obtained using the proposed model 
are validated versus experimental data and compared with numerical results obtained using the standard 
FPV model. 


II. The flamelet—progress-variable models 

For the case of non-premixed combustion of interest here, the basic assumptions of the flamelet model are 
fulfilled for sufficiently large Damkohler number, Da- In fact, when the reaction zone thickness is very thin 
with respect to the Kolmogorov length scale, turbulent structures are unable to penetrate into the reaction 
zone and cannot destroy the laminar flame structure. Effects of turbulence only result in a deformation 
and straining of the flame sheet and locally the flame structure can be described as function of the mixture 
fraction, Z, the scalar dissipation rate, %, and the time. The scalar dissipation rate is a measure of the 
gradient of the mixture fraction representing the molecular fluxes of the species towards the flame and is 
defined as y = 2Dz{^Z)'^, where Dz is the molecular diffusion coefficient of the chemical species. Therefore, 
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the entire flame behavior can be obtained as a combination of solutions of the laminar flamelet equation. In 
the present work we consider a further simplification assuming a steady flamelet behavior, so that chemical 
effects are entirely determined by the value of Z, whereas x describes the effects of the flow on the flame 
structure according to the following steady laminar flamelet equation (SLFE) for the generic variable ip: 

In equation o, p is the density and is the source term related to (p [H], different from zero in the case 
of finite rate chemistry. In particular, in this work, the FPV model proposed by Pierce and Moin [inillH] is 
employed to evaluate all of the thermo-chemical quantities involved in the combustion process. This approach 
is based on the parameterization of the generic thermo-chemical quantity, (p, in terms of the mixture fraction, 
Z, and of the progress parameter. A, instead of y: 


cP = F^{Z,A). 


( 2 ) 


In principle a transport equation for A can be solved, however, it has several unclosed terms that need to 
be modeled m- So that, it is defined a progress variable, C, obtained trough: C = Fc{Z,K). A transport 
equation for C is solved and the flamelet library is so parameterized in terms of Z and C assuming that the 
Fc, see equation ((2), is invertible, 

K = Fc\Z,C). (3) 

Equation ([2]), in terms of Z and C, reads: 

cP = F^iZ,Fc\Z,C)). (4) 

The choice of the progress variable is not unique and some recent works discuss in details this issue proposing 
a procedure for its optimal selection [20U22] . A suitable definition for the progress variable is the sum of the 
mass fraction of the main products |20) ; for hydrogen combustion: 

C = Yh,o- (5) 


Equation ([2]) is taken as the solution of the SLFE ([T|) . The solution variety over y = Xst is called S-curve. A 
key difficulty to integrate the flamelet equation is to know a priori informations about the scalar dissipation 
rate dependence on the mixture fraction, x(^). In fact, flamelet libraries are computed in advance and are 
assumed to be independent of the flow field. Therefore, the dependence of x on Z has to be modeled [23H25] . 
In this work, the functional form of xi^) has been taken from an idealized flow configuration, as proposed by 
Peters [26]; the distribution of the scalar dissipation rate in a counterflow diffusion flame, $(Z), is employed, 
scaled in the following way: 


x{Z) = Xst 


HZ) 

HZst)' 


( 6 ) 


where Xst and Zst are evaluated at the stoichiometric point |26j . From equation (l2|) one can obtain the 
Favre-averaged value of cp and of its variance using the definitions: 


4 > 


F^{Z,A)P{Z,A)dZdA, 


( 7 ) 


A"2 


{F^{Z, A) - (p)^P{Z, A)dZdA, 


where P(Z, A) is the density-weighted PDF, 


P(Z,A) = 


pP{Z,A) 


( 8 ) 


(9) 


P(Z, A) is the joint PDF and p is the Reynolds-averaged density. As usual, (p can be decomposed as: 


= 6 +: 


p^ 

V 


and p = p + p', 


( 10 ) 
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where (f)" and p' are the fluctuations. This function plays a crucial role in the definition of the model, affecting 
both its accuracy and computational costs. Moreover, the choice of such a function is not straightforward 
because of the unknown statistical behavior of the two variables Z and A [H]. The definition of P{Z,A) 
is still an open problem whose solution is being pursued by several researches p7H5D] . The aim of this 
work is to provide a more general model based on the statistically most likely distribution (SMLD) [T7] 
for the joint PDF of Z and A [JT]; this model evaluates the statistical correlation between such variables 
and, using the total energy conservation to evaluate the temperature field, represents a effective method 
for the simulation of compressible reacting flows. Here, the performance of such a combustion model are 
assessed by employing a Reynolds-Averaged Navier-Stokes solver with k-ui turbulence model m to compute 
a hydrogen-air supersonic combustion. 


A. Presumed probability density function model 

In this section, the standard FPV model [TH] (called here Model A) and the joint-PDF SMLD model (called 
here Model B) [3T] are briefly described. 

Since the mixture fraction, Z, and the progress parameter, A, are the two independent variables of the 
combustion model, they form a basis by which one can derive all of the thermo-chemical quantities. When 
a turbulence model is used, the evaluation of the average quantities, see equations © and ([5]), requires the 
PDF to be known or somehow presumed. Such a PDF establishes the statistical correlation between Z and 
A. Employing the Bayes’ theorem, 

P(Z,A)=P(Z)P(A|Z), (11) 

one usually presumes the functional shape of the marginal PDF of Z and of the conditional PDF of A|Z. In 
model A, the basic assumption is the statistical independence between Z and A, so that P{Z, A) = P{Z)P{A). 
Furthermore, the statistical behavior of the mixture fraction is described by a /3 - distribution. In fact, even 
though the definition of P{Z) is still an open question [T^, it has been shown by several authors that 
the mixture fraction behaves like a passive scalar whose statistical distribution can be approximated by a 
P function [32II34] . The two parameter family of the /3-distribution in the interval x G [0,1] is given by: 

/3(x; X, x"^) = x°~^(l - x)'’-^ 

r(a)r(6) 

where F(x) is the Euler function and a and b are two parameters related to x and x"^ 

x(x —x^—x"2) , (1 — x)(x — x^ — x"2) 

a= --, 0 =--. (13) 

^1/2 ^1/2 

Moreover, P(A) is chosen as a Dirac distribution, implying a great simplification in the theoretical framework. 
With these criteria, the Favre-average of a generic thermo-chemical quantity is given by: 


</> 


F^{Z,A)p{Z)5{A - A)dZdA = 


F4Z,A)p{Z)dZ. 


(14) 


Therefore, one has to solve only three additional transport equations (for Z, and A) to evaluate all 
of the thermo-chemical quantities, thus avoiding the expensive solution of one transport equation for each 
chemical species. 

Model Debased on the SMLD approach to model the joint PDF, does not need any assumption about 
the form of P(Z, A). Following such an approach, the probability distribution can be evaluated as a function 
of an arbitrary number of moments of Z and A. It is noteworthy that, even though equation ([T]) is based on 
the assumption that Z and A are independent, one can properly take into account the statistical correlation 
between Z and A employing the SMLD joint-PDF in the evaluation of the effects of turbulence m- 

In this work, the first three moments of the joint probability density function P(x), where x = (Z, A)^, 
are assumed to be known; therefore, the joint-PDF reads m- 


Psml,2{Z,A) = — exp|- Li,i(Z - Z)+ Aii, 2 (A - A) 
Un t L 


/io 

1 

2 L' 


M2,ii(-^ — ^)^ + M2,i2(-^ — ■^)(A — A)-I-^2,2 i(A — A)(Z — Z)-I-/r2,22(A — A)^ |. (15) 
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In the equation above, /ig is a scalar, /Ti is a two - component vector, and |t 2 is a square matrix of rank two: 

Mo = y dxPsML,2i^), (16) 

-Ml,* = y'rf£a,,PsML.2(f)=/3(l;f*,^f)-/3(0;e*,Cf), (17) 

Ski - M2,fen Cn^i = J dxdx^iixi - '^l)PsML,2{x)) = ^(l;Cfe)CfeCi) “ 6:Ml,q (18) 

where z, k, n, and I indicate the vector components; ^i, and are the mean (£)), the fluctuation {xi — Xi) 

and the variance {x'P) of the z-th component of x, respectively; finally, P indicates the beta distribution 

function. _ 

In Model B, one has to solve four additional transport equations (for Z, C, and C"^) to evaluate 
all of the thermo-chemical quantities. 

III. Governing equations 


A. Thermodynamics 


The Peng-Robinson (PR) Equation of State (EoS) [33[3n] is employed in the present work: 

RT a 

Vm-b~ (E^+2Kn6-62)’ 


(19) 


where R is the universal gas constant and Vm is the molar volume. The parameters a{T) and b account 
for the effects of attractive and repulsive forces between molecules, a proper temperature dependence of a 
being essential for predicting vapor pressures correctly [37) . The PR EoS is one of the most frequently used 
cubic equations of state, due to its straightforward implementation and its accuracy [35]. The extension to a 
multi-component mixture, according to the mixing rules proposed by Reid et al. [36] , provides the following 
expressions of the parameters a and b: 


Ns Ns Ns 

a = '^'^XiXjaij(T), b = ^Xth, 

i j i 


( 20 ) 


where is the number of species, Xi is the molar fraction of the species z, and the coefficients aij{T) and 
bi are given in Refs ()381I39[I and are omitted here for brevity. 

To ensure the self-consistency of the model, all of the thermodynamic properties of the flow are calculated 
from the same equation of state. The properties of interest for the present fluid dynamic simulations are the 
specific enthalpy, h, and the constant-pressure specific heat, Cp. Such properties can be obtained from the 
Gibbs energy, G, 

G(T,p) = / p{V^,T,x^)dV^ +pVm -RT + y^x^ [Cf + RT\n{x^)] , (21) 


where the superscript “ref” indicates the low-pressure reference condition [40] corresponding to = 
100kPa, and 14*,„ = i?T/(p"®f): 


h = G-T 



= + pVm - RT + Ki 

p,x 


a-T 


da 

df 


( 22 ) 


and 


Cry - 


dTJ 


= cf - T 


{dpidT)\ 


Ki = 


2V2b 


In 


{dp/dVm)T 
14* + (1 - V2) b 
Vm -b (l -b V^) b 


- - R-T^Ki 


dT-i 


(23) 

(24) 


In equations (I22p and (1331) the partial derivatives of a with respect to the temperature are evaluated analyt¬ 
ically, see Ref (1301) for their expressions. 
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B. Transport properties 

Two advanced models have been implemented in this work to evaluate correctly the viscosity m and the 
thermal conductivity [42] , which strongly depend on the pressure at the thermodynamic conditions of interest 
here. 

Cho and Chung [41] , starting from a base formulation valid at low pressures, developed an expression for 
the viscosity which is valid at high (supercritical) pressure and low (transcritical) temperature: 


M M 


36.344 




2/3 


(25) 


see Ref ()36[l for the meaning of the symbols and for further details. 

Ely and Hanley [42] used Eucken’s splitting of the thermal conductivity into contributions from the 
interchanges of both translational and internal energies [36] . together with a corresponding state method 
using methane as the reference component to estimate the thermal conductivity of non-polar fluids over a 
wide range of densities and temperatures: 


^ + '^ Xi Xj I 

i j 

see again Ref (l36l) for the meaning of the symbols and for further details. 


(26) 


C. Flow equations and numerical solution procedure 

The numerical method developed in Ref (ITHl) has been employed to solve the steady-state RANS equations 
with k-uj turbulence closure. For an axisymmetric multi-component reacting compressible flow the system 
of the governing equations can be written as: 

dtQ + d,{E-E,)+dy{F-F,) = S, (27) 


where t is the time variable; x and y are the axial and the radial coordinate, respectively; k and w are the 
turbulence kinetic energy and its specific dissipation rate; Q={'p,'pUx,'pUy,'pH — pt,'pk, pa;,pi?„) is the vector 
of the conserved variables; p, (ux, Uy), H indicate the Reynolds-averaged value of density, the Favre-averaged 
values of velocity components and specific total enthalpy given hy H = h+^{ii^ + Uy) + ^k with h accounting 
for the species enthalpy per unit mass, respectively; is a generic set of conserved variables related to the 
combustion model. In this framework, is the set of independent variables of the flamelet model, namely, 
Z, C, (see the following subsection); E, F, and Fy are the inviscid and viscous flux vectors [15] . 
respectively; S is the vector of the source terms. 

The heat flux in the total energy equation, namely q = —pCpEx'^T + pVnYnhn, is composed by 

two terms since the Dufour effect is neglected, Dt is the thermal diffusivity and Cp the specihc heat at 
constant pressure. The mass diffusion term is treated with the Fick’s law considering Vn = 
assigning a mixture diffusivity, Dn,mix, to each species. A cell-centered finite volume space discretization 
is used on a multi-block structured mesh. The convective and viscous terms are discretized by the third- 
order-accurate Steger and Warming [44] flux-vector-splitting scheme and by second-order-accurate central 
differences, respectively. An implicit time marching procedure is used with a factorization based on the 
diagonalization procedure of Pulliamm and Chaussee [45] . employing a scalar alternating direction implicit 
(ADI) solution procedure [46]. Steady flows are considered and the ADI scheme is iterated in the pseudo¬ 
time until a residual drop of at least five orders of magnitude for all of the conservation-law equations is 
achieved. Characteristic boundary conditions for the flow variables are imposed at inflow and outflow points, 
whereas no slip and adiabatic conditions are imposed at walls; fc, uj, and are assigned at inflow points, 
whereas they are linearly extrapolated at outflow points. At solid walls, k is set to zero and uj is evaluated 
as proposed by Menter and Rumsey [47]: 


w = 60 


0.09 


(28) 


where is the distance of the first cell center from the wall; the homogeneous Neumann boundary condition 

is used for (non-catalytic wall). Finally, symmetry conditions are imposed at the axis. 
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Table 1. Conditions for RCM-3 2001 test-case. 


rh (kg/s) T (K) p (kg/m^) u (m/s) Tui„ (-) At, in, mm 
LO 2 0.1 85 1177.8 4.35 5% 4 

H 2 0.07 287 5.51 236 5% 4 


D. Turbulent FPV transport equations 

For the case of turbulent flames, the solution of the SLFE, namely equation ([2]), is expressed in terms of the 
Favre averages of Z and C and of their variances. Using Model A, one can tabulate all chemical quantities 
in terms of Z, Z'"^ and C, since the model is independent of C'"^. On the other hand, Model B expresses 
(j) also in terms of C'"^ and therefore an additional transport equation needs to be solved. In this case, the 
transport equations for the combustion model (included in equation (1271) 1 are written as: 


dt{pZ) + V • {puZ) 
5t(p^) + V • (p5^) 
dtipC) + V • (p^C) 
dt{p^^) + V • (p5c^) 


V • 

{p + U|)pVZ 



(29) 

V- 


-pX + 2pU|(VZ)2, 

(30) 

V • 

{D + 

+ pwc. 

(31) 

V- 


- pxc + 2-pD^pVC f + 2pC^, 

(32) 


where xc is modeled in terms of Z"'^ and C"^ [28], namely xc = ^ i® i'ii® diffusion coefficient for all 

of the species, given as D = vjPr assuming a unity Lewis number; v and Pr are the kinematic viscosity 
and the Prandtl number, respectively; D*~ = = D*~ = 7?'^ = vtjSct are the turbulent mass diffusion 

coefficients, Sct being the turbulent Schmidt number; finally, (oc is the source term for the progress variable 
precomputed and tabulated in the flamelet library. The gradient transport assumption for turbulent fluxes 
is used and the mean scalar dissipation rate, x s-nd Xc, appear as a sink term in equations (1501) and (I32L 
respectively. 

At every iteration, the values of the flamelet variables are updated using equations (l29]) - (l32l) and the Favre- 
averaged thermo-chemical quantities are computed, using equation (ED. Such solutions provide the mean 
mass fractions which are used to evaluate all of the transport properties of the fluid, namely the molecular 
viscosity, the thermal conductivity and the species diffusion coefficients. 


IV. MASCOTTE VOS: H2/LO2 supercritical combustion 

The MASCOTTE cryogenic combustion test facility was developed by ONERA to study fundamental 
processes in cryogenic combustion involving propellants such as liquid oxygen and gaseous hydrogen. In 
particular, the test-case RCM-3 2001 [IS], dealing with the super-critical H 2 /LO 2 combustion problem, has 
been chosen here as a suitable test for investigating: 1) the influence of the real-gas model on the simulated 
combustion phenomenon; 2) the importance of using an accurate chemical kinetic scheme within the flamelet 
model; 3) the influence of the presumed probability function model on the flame structure. 

The injector consists of an inner diverging duct for the oxygen with inlet diameter of 3.6 mm and outlet 
diameter of 5 mm. Hydrogen is injected coaxially, through an annular duct with inner and outer diameters 
equal to 5.6 mm and 10 mm, respectively. The injector length is equal to 50 mm so that a fully developed 
turbulent flow is obtained at the injector exit. The combustion chamber has a 50 x 50 mm^ square section, 
the edge length being 50 mm. In order to perform an axisymmetric simulation, the combustion chamber is 
modeled as a cylinder with a radius of 28.21 mm such as to preserve the chamber section area. 

The chamber pressure was held at 6 MPa, higher than the critical pressure for both oxygen and hydrogen 
(5.043 MPa and 1.313 MPa, respectively). Liquid oxygen is injected at the temperature of 85 K, correspond¬ 
ing to a density of 1177.8 kg/m^, whereas hydrogen is injected at a temperature of 287 K. A computational 
grid with about 18,000 cells distributed among 24 structured blocks has been used. A local view of the grid 
is given in figure |T| Table |T] provides the inlet boundary conditions. 
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Figure 1. Reference geometry for Test RCM-3 2001 and a detail of the computational grid for the injectors. 


For this test case, a limited experimental dataset for comparison is available. Quantitative experimental 
profiles of temperature are available at three axial locations {x =0.015, 0.05 and 0.1m), with a quite 
poor resolution (four points at r = 0.004, 0.008, 0.012 and 0.016m). On the other hand, any quantitative 
experimental data are available for chemical species. Additionally, a LIF image data m of the hydroxyl 
radical is used for qualitative comparison. 

Steady flamelet calculations have been performed with the chemical kinetic scheme provided by Li et 
al. |51) . from now on called Li-scheme, who developed a scheme consisting of 25 elementary reactions and 9 
species, which was proven accurate for operating pressures up to 9 MPa. In the present study, the chemistry 
library is discretized with 125 uniformly distributed points in the Z and C directions, and 25 points in the 
Z''^ and C""^ direction. 


A. Gas model effects 


The system of the flamelet equations is solved here employing either the ideal-gas EoS, with state-of-the-art 
models for the temperature dependence of the thermodynamic and transport properties |52] , or the real-gas 
EoS with the advanced transport properties described in sections El and IS 

The influence of the real-gas model on the flame structure, and consequently on the flamelet calculation, 
is quite remarkable, as just pointed out by some authors [T6l[53]. The application of the real-gas model 
reduces the maximum temperature, due to a lower heat release, and consequently predicts lower peaks for 
almost all of the species concentrations. For an H 2 /L 02 combustion at a 6 MPa pressure, the two S-shaped 
curves computed by using either the ideal- or the real-gas models are provided in figure[2] Due to the higher 
dissipation rate, the ideal-gas model provides a higher hydrogen consumption rate along the entire stable 
upper branch, resulting in higher flame temperatures, as shown in figure [ ^a)[ Accordingly, a higher heat 
release is predicted by the ideal-gas flamelet with respect to the real-gas one, with the peak shifted towards 
the oxidizer side, see figure [ ^b)[ Therefore, the lower fuel consumption rate for the real-gas flamelet leads 
to lower species concentrations, as shown in figures [ ^a)|(b)' the mole fractions for the final combustion 
products and radical species are lower for the real-gas flamelet for all species except HO 2 . 

These sets of flamelet have been integrated accordingly to both Model A and Model B approaches, but 
only the first will be used in this analysis on real-gas effects. 

Firstly, a simulation has been performed using the ideal-gas model for both the equation of state and 
the flamelet model. The ideal-gas equation of state fails in the prediction of the density of oxygen at the 
injection temperature of 85 K and a chamber pressure of 6 MPa: it provides a value of 271.7 kg/m^ versus 
an experimental value [54] of 1177.8 kg/m^. The main issue of such model is that the computed oxygen 
mass-flow rate drops to 0.023 kg/s. 
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Figure 2. S-shaped curve obtained with ideal- (solid line) and real-gas (dashed line) models. 



Figure 3^ _ Ideal- (solid lines) and real-gas (dashed 

density; |(b)| heat release. 



Mixture fraction, Z 


(b) 

solutions at Xst = 10® s“^. |(a)| temperature and 




Figure 4. Ideal-gas (solid line) and real-gas (dashed line) solutions at Xst = 10® s~^. |(a)] final products: H 2 O 
(plain line), HO 2 ■ 100 (open circles), H 2 O 2 ■ 100 (solid diamonds); [(b)] radical species: OH (plain line), O (open 
circles), H (solid diamonds). 


Figure [S] provides a comparison between numerical results and experimental data for OH mass fractions 
(left), and a comparison between numerical results for temperature (right). As expected, the ideal-gas model 
(a) does not appear adequate for the present test-case, providing a very short flame, essentially due to the 
reduced oxygen mass-flow-rate, and consequently a different oxygen/fuel ratio. 

Then, the influence of the gas-model onto the derivation of the flamelet library has been investigated. 
The influence of the gas-model on the temperature field is quite remarkable, as shown in figure [SJ which 
provides the results obtained using the ideal- and real-gas models in the flamelet approach, indicated by 
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Figure 5. OH mass fractions (left) and temperature (right): (top-left) Abel-transformed-emission image, (a) 
Ideal-gas EoS — Ideal-gas flamelet model, (b) Real-gas EoS — Ideal-gas flamelet model, (c) Real-gas EoS — 
Real-gas flamelet model. 
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Figure 6. Radial temperature distributions at several axial locations: Ideal-gas EoS — Ideal-gas flamelet (solid 
line), Real-gas EoS — Ideal-gas flamelet (dashed line), Real-gas EoS — Real-gas flamelet (dashed-dotted line) 
models 
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Figure 7. Ignition delay of H2/02/Ar mixtures in shock tubes. Li-scheme (dashed-dotted-dotted line), Warnatz- 
scheme (solid line), and Jachimowski-scheme (dashed line). 


label (b) and (c), respectively. It appears that, due to the higher dissipation rate, the ideal-gas flamelet 
model provides a higher hydrogen consumption rate and, therefore, a shorter flame, which seems closer to 
the experimental result that provides a closed flame at about x = 0.09m. However, the ideal-gas flamelet 
model exhibits a higher spreading angle, in particular close to the injector, where the flame should be a 
thinner layer confined around the oxygen core. 

A more detailed comparison is given in figure IHl showing the radial temperature distributions at six 
positions along the axis. Close to the injector {x = 0.015m), the real-gas flamelet model results are in 
reasonable agreement with the experimental data except the hot flame zone close to the axis, where CARS 
data are not sufficiently resolved to capture the flame layer. On the other hand, the ideal-gas flamelet 
model provides a wider flame due to the higher spreading angle observed. Slightly apart from the injector 
(x = 0.05m) both the ideal- and real-gas flamelet model overestimate the temperature in the hydrogen-rich 
flame zone (r = 0.012 and 0.016m), but again, any indication on the position of the flame front (somewhere 
between r = 0.004 and 0.012m) is provided by experimental data which exhibits a flat profile of about lOOOK. 
Finally, radial distributions confirm that the ideal-gas flamelet model provides a faster combustion, leading 
to low values of the temperature gradient, characteristic of the post-combustion diffusion process, already 
aX, X = 0.1 TO, where experimental profile indicates that the hot flame zone could be placed in the proximity 
of the jet core. 

B. Kinetic scheme effects 

Here, in order to assess the role of the basic chemical kinetic scheme employed for the definition of the 
steady flamelet equations, three alternative schemes have been considered. Two of them are simplified 
models proposed by Jachimowski |55] and Marinov [56) . from now on called Jachimowski- and Marinov- 
scheme, respectively; the former is a reduced scheme with 7 reactions and 7 species, whereas the latter is a 
global scheme (3 species). The third one is an additional detailed schemes provided by Warnatz [57] (from 
now on called Warnatz-scheme). 

The performance of these mechanisms are at first evaluated in a shock tube laminar test for which 
experimental data, namely ignition delay times, are available [ 55 ]. The AURORA code [ 55 ], was used to 
simulate experimental conditions in the shock tube. The flow is assumed zero-dimensional, laminar and 
well mixed. The operating pressure is p = 6.5 MPa and the initial composition is set to: 1(32 = 0.01, 
Yhs = 0.02, and Yat = 0.97, whereas the initial temperature ranges from 1000 to 2000 K. The ignition 
delay time is defined as the time when [OH]/t reaches its maximum, where [OH] = pToh/AIoh indicates 
the molar concentration of the species OH. As a consequence of this definition, the ignition delay time 
cannot be calculated with the Marinov-scheme, being the OH species concentration not calculated. In 
figure [3 the ignition delay times of H2/02/Ar mixtures are reported for the Li-, Warnatz- and Jachimowski- 
scheme. It appears that (1) the reduced Jachimowski-scheme provides very short ignition delay at low initial 
temperatures, diverging from the experimental at about 1400 K; (2) the detailed schemes are able to capture 
the rapid increase of the ignition delay at lower temperature, that as indicated by Refs ([SKSO]), is due to 
the increased chain-termination mechanisms. 
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Then, the considered schemes have been used for the generation of flamelet libraries and 2D axisymmetric 
CFD computation have been performed. All the results have been obtained using the real-gas flamelet model 
and the turbulent combustion Model A. 


Kinetic scheme effects 
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Figure 8. OH mass fractions (left) and temperature (right): (top-left) Abel-transformed-emission image, (a) 
Marinov-scheme, (b) Jachimowski-scheme (c) Warnatz-scheme, (d) Li-scheme. 


Figure [5] shows OH mass fractions and temperature contours obtained using the Marinov-, Jachimowski-, 
and Warnatz-schemes, respectively, the results computed using the Li-scheme being given for comparison, 
too. The results obtained with the simplified scheme are unacceptable, whereas the two detailed schemes 
qualitatively reproduce the experimental distributions. A more detailed comparison is provided in figure IHl 
which shows the temperature distributions at several axial locations obtained using the four kinetic schemes. 
The results obtained using the two simplified schemes indicate again that the flame is very short with respect 
to the experimental data as well as to the numerical results obtained using either detailed scheme, which 
appear to be more accurate. 

In fact, both simplified schemes predict a barely appreciable ignition delay, the combustion developing as 
soon as the reactants come into contact. Therefore, the combustion takes place very close to the injector and 
the reactants cannot be transported by the flow further downstream. Such results confirm that, especially for 
high-pressure combustion processes, a detailed kinetic scheme is warranted, ft is noteworthy that different 
temperature distributions are provided by the Li- and Warnatz-scheme, the latter predicting a slightly shorter 
flame with a smaller spreading angle and a thinner reaction zone. 

C. Turbulent combustion model effects 

The results shown in the previous sections have been obtained using the flamelet Model A, which based on 
the original flamelet-progress-variable model developed by Pierce et al. [18]. Here, in order to asses the role 
of the presumed probability function on high-pressure combustion, the flamelet Model B has been used to 
simulate the MASCOTTE test case. All the results of this section have been obtained using the real-gas 
flamelet model and the Li-scheme. 

Figure [To] shows OH mass fractions and temperature contours obtained using the Model A and Model B. 

The results obtained by Model B seems qualitatively in better agreement with experimental data with 
respect to Model A. At first. Model B shows a shorter flame length with respect to Model A and this seems 
to be in accordance with experimental observations. Moreover, the characteristic bump of the flame front 
due to the entrainment of cold fluid from the corner vortex is moved backward, likewise to the experimental 
observation, indicating that the injection and mixing processes combined with combustion are better cap- 
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Figure 9. Radial temperature distributions at several axial locations: Marinov-scheme (solid line), 
Jachimowski-scheme (dashed line), Warnatz-scheme (dash-dotted line), Li-scheme (dotted line) models. 
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Figure 10. OH mass fractions (left) and temperature (right): (top-left) Abel-transformed-emission image, (a) 
Model A, (b) Model B. 
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Figure 11. Radial temperature distributions at several axial locations: Model A (dashed line), Model B (solid 
line). 


tured. Last but not least, the flame thickness provided by Model B grows more slowly in the first part of 
the flame close to the injector, where the flame should be a very thin hot layer around the cold L02 core, 
whereas it became larger moving slightly forward. 

A more detailed comparison is provided in figure [5] which shows the temperature distributions. The 
results indicate again that the Model B flame is very thin close to the injectors {x = 0.015m) and computed 
temperature is in very good agreement with experimental data. In the section at a: = 0.1m, Model B 
shows a better agreement with the experimental data, the temperature is underestimate but the behavior 
is well predicted. It seems, in fact, that the combustion process is too fast in the last part of the flame. In 
conclusion, seems that the more general framework given by Model B can actually improve the prediction 
capabilities of the combustion model joined with the real gas equation of state. The greatest differences are 
observed considering the two flame shapes and the initial part, near to the injector x <0.03 m, due to the 
better evaluation of the corner vortex. 

D. Conclusions 

This paper provides a numerical method based upon RANS equation for the simulation of high-pressure 
conditions. The turbulent combustion coupling has been modeled by implementing flamelet-progress-variable 
models, with the H 2 /L 02 combustion kinetics provided by four kinetic schemes (both reduced and detailed), 
and thermodynamics by the Peng Robinson real-gas equation of state. Moreover, a general framework for 
the evaluation of the most probable joint distribution of the mixture fraction and the progress variable has 
been developed and used to generate a flamelet lookup-table. 

The MASCOTTE VOS test case has been computed, involving the supercritical combustion of H 2 /L 02 
and the effects of the modeling on the results has been analyzed. In particular, the test: 

• allowed for a detailed assessment of the real-gas effects, on both the main flow equations and flamelet 
calculations; 

• showed the importance of using a detailed kinetic scheme; 
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• proposed a more general framework for the derivation of the density probability function for FPV 
models that could improve actual predictions. 
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